#region Preset using directives

using System;
using System.Collections.Generic;
using System.Drawing;
using System.IO;
using System.Linq;

using VirtualLab.Programming;
using VirtualLabAPI.Core.BasicFunctions;
using VirtualLabAPI.Core.Common;
using VirtualLabAPI.Core.DataVisualization;
using VirtualLabAPI.Core.FieldRepresentations;
using VirtualLabAPI.Core.Functions;
using VirtualLabAPI.Core.GeometryDescription;
using VirtualLabAPI.Core.LightPath;
using VirtualLabAPI.Core.Materials;
using VirtualLabAPI.Core.Numerics;
using VirtualLabAPI.Core.Numerics.Region2D;
using VirtualLabAPI.Core.OpticalSystems;
using VirtualLabAPI.Core.Propagation;
using VirtualLabAPI.Core.SupportFunctions;
using Wyrowski.Infrastructure.Support;

#endregion

namespace OwnCode {

	/****************************************************************************
	Title: "Design Meta Structure from Target Function"
	Author: ChHe (LightTrans)
	Description: This module designes a pillar distribution for a meta surface by using a desired phase function and
	the look-up table of the complex characteristic of the meta atoms (e.g. pillars) over e.g. their diameter.
	Outputs are the distribution of pillars, the complex response of the designed element based on the chosen pillars
	and the deviation between th complex response and the initially desired phase.

	Version (Date): 1.0 (2023-06-20)
	****************************************************************************/
	public class VLModule {
		//define amplitude threshold which shall be used to select the phase values
		//value is given as multiple of the maximum amplitude
		//if the amplitude value is smaller than the threshold the set parameter is not used for phase design
		private double amplitudeThreshold = 0.9;

		public void Run() {
			#region get data
			//ask for the target phase function
			DataArray2D targetFunction = (DataArray2D)VL_GUI.SelectOpenDataArray2D("Select Target Phase Function");
			//error handling (no target function selected)
			if (targetFunction == null) {
				Globals.DataDisplay.LogError("No target function selected.");
				return;
			}

			//error handling (only one subset)
			if (targetFunction.Data.Length > 1) {
				Globals.DataDisplay.LogError("Selected target function contains more than one subset!");
				return;
			}

			//error handling (is complex valued)
			if (targetFunction.IsComplex) {
				Globals.DataDisplay.LogError("Selected target function is complex valued! Real-valued phase information (in 2π modulo) expected!");
				return;
			}

			//ask for the lookup table to perform the pillar design
			DataArray1D lutComplexResponseVsDiameter = (DataArray1D)VL_GUI.SelectOpenDataArray1D("Select Lookup Table Phase vs. Diameter");

			//error handling (no LUT selected)
			if (lutComplexResponseVsDiameter == null) {
				Globals.DataDisplay.LogError("No lookup table selected");
				return;
			}

			//error handling (only one subset)
			if (lutComplexResponseVsDiameter.Data.Length > 1) {
				Globals.DataDisplay.LogError("Selected lookup table contains more than one subset!");
				return;
			}

			//error handling (is complex valued)
			if (!lutComplexResponseVsDiameter.IsComplex) {
				Globals.DataDisplay.LogError("Selected lookup table is not complex valued!");
				return;
			}
			#endregion

			//generate matrix to store the pillar diameters
			CFieldDerivative2DReal cfPillarDiameters = new CFieldDerivative2DReal(new Vector(targetFunction.NoOfDataPoints_X,
																							 targetFunction.NoOfDataPoints_Y));
			//generate matrix to store the response function evaluated from the lookup table
			ComplexField cfResponseOfMetaElement = new ComplexField(new Vector(targetFunction.NoOfDataPoints_X,
																			   targetFunction.NoOfDataPoints_Y));
			//generate matrix to store the phase residual (difference from target phase and resulting designed phase)
			CFieldDerivative2DReal cfPhaseResidual = new CFieldDerivative2DReal(new Vector(targetFunction.NoOfDataPoints_X,
																						   targetFunction.NoOfDataPoints_Y));



			//detect minimum and maximum amplitude within LUT
			double minAmplitude, maxAmplitude;
			lutComplexResponseVsDiameter.Data[0].Normalize();
			ComplexFieldEvaluation.GetMinMaxValues(lutComplexResponseVsDiameter.Data[0],
												   FieldQuantity.Amplitude,
												   out minAmplitude,
												   out maxAmplitude);

			//loop over all data points of the target function
			for (int runY = 0; runY < targetFunction.NoOfDataPoints_Y; runY++) {
				for (int runX = 0; runX < targetFunction.NoOfDataPoints_X; runX++) {
					//extract phase value from target function
					double targetPhaseOnCurrentCell = targetFunction.Data[0][runX, runY].Re;
					//define index within LUT that shall be used to realize current target phase
					int foundIndex = -1;

					#region select pillar diameter
					//define min phase value for comparison
					double minPhaseDiff = double.MaxValue;

					//loop over all LUT entries
					for (int runLUTEntries = 0; runLUTEntries < lutComplexResponseVsDiameter.NoOfDataPoints; runLUTEntries++) {
						//get current value
						Complex currValueLUT = lutComplexResponseVsDiameter.Data[0][runLUTEntries];

						if (currValueLUT.Abs() > this.amplitudeThreshold * maxAmplitude) {
							//get phase from LUT entry
							double phaseMetaAtom = currValueLUT.Arg();
							//calculate phase different
							double phaseDiff = Math.Abs(targetPhaseOnCurrentCell - phaseMetaAtom);

							//check whether phase of current meta atom is better than the already found
							if (phaseDiff < minPhaseDiff) {
								//store index
								foundIndex = runLUTEntries;
								//store minimum found phase difference
								minPhaseDiff = phaseDiff;
							}
						}
					}
					#endregion

					if (foundIndex < 0) {
						throw new ArgumentException("No value found for target phase.");
					}

					//extract information for found LUT entry
					double pillarDiameterToUse = lutComplexResponseVsDiameter.CoordinateOfDataPointIndex(foundIndex);
					Complex selectedMetaResponseValue = lutComplexResponseVsDiameter.Data[0][foundIndex];

					//store value of pillar
					cfPillarDiameters[runX, runY] = pillarDiameterToUse;
					//store response value of meta atom
					cfResponseOfMetaElement[runX, runY] = selectedMetaResponseValue;
					//calculate and store phase difference
					cfPhaseResidual[runX, runY] = targetPhaseOnCurrentCell - selectedMetaResponseValue.Arg();
				}
			}

			//generate ouput for pillar diameter distribution
			DataArray2D daPillarDistribution = new DataArray2D(new ComplexFieldArray(new ComplexField[] { cfPillarDiameters }),
															   new PhysicalProperty[] { PhysicalProperty.Length },
															   new string[] { "Pillar Diameter" },
															   targetFunction.SamplingDistance_X,
															   targetFunction.CoordinateOfFirstDataPoint_X,
															   targetFunction.PhysicalPropertyOf_X_Coordinates,
															   targetFunction.CommentOfCoordinates_X,
															   targetFunction.SamplingDistance_Y,
															   targetFunction.CoordinateOfFirstDataPoint_Y,
															   targetFunction.PhysicalPropertyOf_Y_Coordinates,
															   targetFunction.CommentOfCoordinates_Y);
			//display document
			Globals.DataDisplay.ShowDocument(daPillarDistribution, "Pillar Distribution");


			//generate ouput for response function of the meta element
			DataArray2D daResponseFunction = new DataArray2D(new ComplexFieldArray(new ComplexField[] { cfResponseOfMetaElement }),
															 new PhysicalProperty[] { PhysicalProperty.ElectricalField },
															 new string[] { "Complex Response" },
															 targetFunction.SamplingDistance_X,
															 targetFunction.CoordinateOfFirstDataPoint_X,
															 targetFunction.PhysicalPropertyOf_X_Coordinates,
															 targetFunction.CommentOfCoordinates_X,
															 targetFunction.SamplingDistance_Y,
															 targetFunction.CoordinateOfFirstDataPoint_Y,
															 targetFunction.PhysicalPropertyOf_Y_Coordinates,
															 targetFunction.CommentOfCoordinates_Y);
			//display document
			Globals.DataDisplay.ShowDocument(daResponseFunction, "Response Function of Meta Element");


			//generate ouput for phase residual of meta element design
			DataArray2D daPhaseResidualDesign = new DataArray2D(new ComplexFieldArray(new ComplexField[] { cfPhaseResidual }),
																new PhysicalProperty[] { PhysicalProperty.AngleRad },
																new string[] { "Phase Residual" },
																targetFunction.SamplingDistance_X,
																targetFunction.CoordinateOfFirstDataPoint_X,
																targetFunction.PhysicalPropertyOf_X_Coordinates,
																targetFunction.CommentOfCoordinates_X,
																targetFunction.SamplingDistance_Y,
																targetFunction.CoordinateOfFirstDataPoint_Y,
																targetFunction.PhysicalPropertyOf_Y_Coordinates,
																targetFunction.CommentOfCoordinates_Y);
			//display document
			Globals.DataDisplay.ShowDocument(daPhaseResidualDesign, "Phase Residual after Meta Element Design");
		}
	}
}
